function FVlinearadvectionFOU1D
% File name: FVlinearadvectionFOU1D.m
% Description: Solves the single 1D PDE, dU/dt + div(H) = 0,
% using the FOU finite volume scheme on a uniform Cartesian mesh which is
% a single row of cells in the x direction with sides dx and dy.
% This equation corresponds to the finite volume form:
%
% doubleintegral(dU/dt dA) + lineintegral(H.n ds) = 0.
%
% In this case the equation is the linear advection equation so:
% velocity: v = vx i + vy j where vx is constant and vy = 0.
% flux density: H = vU = vx U i
%
% There are NI cells in the i (x) direction and 1 cell in the j (y)
% direction. Since vy = 0 there is flow only through the left and right
% sides of each cell. Cell areas and side vectors are the same for each cell.
%
% The program is written in a structured way using subfunctions so that it
% can be modified easily to the pseudo-1D and 2D cases on non-Cartesian
% structured meshes.
%
% Initial conditions: Gaussian along 1D line of cell centres.
%
% Boundary conditions: Left is Dirichlet, u=0.
%
% subfunctions: freadmesh, fcellarea, fsidevectors, finitialu,
% finterfaceflux, fghostcells, fcalcdt, fplotresults,
% fdrawmesh.
%
% Verification: Looks correct (for vx>0) when comparing graphs of
% numerical and exact solutions.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
clc; clf; clear all;
runtime=2.0; % runtime in seconds
t=0;         % current time
timelevel=0; % current time level
vx=20;       % velocity component in x direction [m/s]
v=[vx,0];    % velocity vector = vx i + 0 j
[x,y,xcen,solid]=freadmesh; % gets coordinates of the 1D computational mesh
% and solid flags (excluding ghost cells)
% Note that there are no solid interior cells in 1D.
disp('mesh read in')
dy=y(1,2)-y(1,1);
[m,n]=size(x);
NI=m-1; % number of computational cells in i direction.
A=fcellarea(x,y); % computes constant cell area
disp('calculated cell area')
S=fsidevectors(x,y); % compute and store cell side vectors
disp('calculated cell side vectors')
u0=finitialu(xcen,t,vx); % Put initial cell centre values of u in a NIx1 array
uinitial=u0; % store initial profile for plotting
disp('inserted initial conditions')
% Append extra cells around arrays to store any ghost values.
u0=finsertghostcells(u0); % u0 is now (NI+2)x1
u1=zeros(size(u0)); % (NI+2)x1 array for u values at next time level
%%
disp('time marching starts')
while(t<runtime)
timelevel=timelevel+1
u0=fbcs(u0); % Implement boundary conditions using ghost cells.
% u0 is (NI+2)x1 and each computational cell is
% indexed i=2 to NI+1.
dt=fcalcdt(A,S,v); % Finds stable time step for each iteration.
t=t+dt % update time
for i=2:NI+1
IH=finterfaceflux(v,u0,i); % gets the left and right interface fluxes for cell i
IHright=[IH(1,1),IH(1,2)]; % interface flux vector right side
IHleft=[IH(2,1),IH(2,2)]; % interface flux vector for left side
%
sright=[S(1,1),S(1,2)]; % side vector for right side
sleft=[S(2,1),S(2,2)]; % side vector for left side
% FV scheme
% compute total flux out of cell i
totalfluxout=dot(IHright,sright)+dot(IHleft,sleft);
% totalfluxout=vx*u0(i,1)*dy+vx*u0(i-1,1)*(-dy);
u1(i,1)=u0(i,1)-(dt/A)*totalfluxout;
end % of i loop
%
u0=u1; % update u values
end % of while loop
disp('1D FOU scheme ended at time t')
t
%fdrawmesh(x,y,solid)
%fdisplay(u0); % print results for a small mesh as a matrix
uexact=finitialu(xcen,t,vx); % exact solution
fplotresults(xcen,u0,uinitial,uexact) % plot results
disp('finished main program, see figures')
end % FVlinearadvectionFOU1D
%%%%%%%%%%%%%%%%%%%%%%%%%%% subfunctions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,y,xcen,solid]=freadmesh
% Verification: not verified
% The mesh is structured and has NI cells in the i (x) direction and
% 1 cell in the j (y) direction.
% The x and y coordinates of the lower left hand corner of cell (i,j) are
% held in arrays x and y respectively which are both (NI+1)x2. In
% this way the 4 vertices of cell (i,j) are (x(i),y(j)), (x(i+1),y(j)),
% (x(i+1),y(j+1)) and (x(i),y(j+1)). solid is an NI by 1 array which
% flags solid cells. If cell (i,j) is solid then solid(i,j)=1 otherwise
% solid(i,j)=0.
%
NI=100;
x=zeros(NI+1,2); % allocate correct sized array
y=zeros(NI+1,2); % allocate correct sized array
xcen=zeros(NI,1); % allocate correct sized array
solid=zeros(NI,1); % allocate correct sized array
dx=2; % cell length in x direction
dy=1; % cell length in y direction (arbitrary)
for i=1:NI+1
	for j=1:2
		x(i,j)=(i-1)*dx;
		y(i,j)=(j-1)*dy;
	end % of j loop
end % of i loop
% find cell centres by averaging coordinates of vertices
% (this works for a general structured mesh)
for i=1:NI
xcen(i,1)=(x(i,1)+x(i+1,1)+x(i+1,2)+x(i,2))/4;
end % of i loop
end % freadmesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function A=fcellarea(x,y)
% Verification: verified.
% In a uniform Cartesian mesh cell area = dx dy
dx=x(2,1)-x(1,1);
dy=y(1,2)-y(1,1);
A=dx*dy
end % fcellarea
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function S=fsidevectors(x,y)
% Verification: not verified
% For each cell:
% right side vector = dy i + 0 j
% left side vector = -dy i + 0 j
%
S=zeros(2,2);
dx=x(2,1)-x(1,1);
dy=y(1,2)-y(1,1);
% Right side vectors
S(1,1)=dy;
S(1,2)=0;
% Left side vectors
S(2,1)=-dy;
S(2,2)=0;
end % fsidevectors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function u=finitialu(xcen,t,vx)
% Verification: verified
% Inserts initial or exact u values at time t.
[NI,NJ]=size(xcen);
u=zeros(NI,1);
% Initial 1D Gaussian function based on the centre, x=xc, of the
% computational domain.
xmax=max(max(xcen)); % max x value in mesh
xmin=min(min(xcen)); % min x value in mesh
xc=(xmax+xmin)/2; % approx x coord of mesh centre
%
cutoff=(xmax-xc)/2; % cut off radius for Gaussian
for i=1:NI
x=xcen(i,1)-vx*t;
d=abs((x-xc)); % distance from centre
if (d<cutoff)
u(i,1)=exp(-0.01*d^2); % Gaussian
else
u(i,1)=0;
end
end % of i loop
end % finitialu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function IH=finterfaceflux(v,u0,i)
% Verification: not verified
% Calculates right and left interface fluxes for each cell.
% Flux H depends on U, i.e. H = H(U).
% For the 1D linear advection equation, H(U) = vx U i + 0 j.
%
% There is considerable choice for interface flux estimation: different
% choices give different FV schemes.
%
%% The following FOU scheme simply takes:
% IHright=H(u0(i,1)) = v u0(i,1)
% IHleft=H(u0(i-1,1)) = v u0(i-1,1)
% This works (is an upwind scheme) for vx > 0.
%
IH=zeros(2,2); % array to store left and right interface flux vectors
%
vx=v(1);
vy=v(2);
% Right
IH(1,1)=v(1)*u0(i,1); % i component of right interface flux
IH(1,2)=0; % j component of right interface flux
% Left
IH(2,1)=vx*u0(i-1,1); % i component of left interface flux
IH(2,2)=0; % j component of left interface flux
end % fIflux
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function outarray=finsertghostcells(inarray)
% Verification: not verified
% Assumes that ghost cells are needed at either end of the domain.
% array is embedded into (m+2)x1 array of zeros. Hence computational
% indices go from i=2 to m+1. In the FOU scheme a ghost cell
% is needed only at the left end.
[m,n]=size(inarray);
dummy=zeros(m+2,1);
dummy(2:m+1,1)=inarray;
outarray=dummy;
end % finsertghostcells
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function outarray=fremoveghostcells(inarray)
% Verification: not verified
% % Removes ghost cells so that a mx1 array becomes (m-2)x1
% [m,n]=size(inarray);
% outarray=inarray(2:m-1,1);
% end % fremoveghostcells
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function u=fbcs(u)
% Verification: not verified
% Implements boundary conditions using ghost cells index by i=1, i=NI+2.
[m,n]=size(u);
NI=m-2;
% Dirichlet: u=0 at each end.
% Don't need to do anything as u values in ghost cells are already zero.
end % bcs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dt=fcalcdt(A,S,v)
% Verification: not verified
% Finds allowable time step dt using heuristic formula.
F=0.9; % tuning factor (F=1 gives exact results)
rightside=[S(1,1),S(1,2)];
dt=A/abs(dot(v,rightside)); % heuristic formula
dt=F*dt;
end % fcalcdt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fdrawmesh(x,y,solid)
% Verification: not verified
% Description: Draws a structured 2D finite volume mesh and fills in any
% solid cells.
% Date structures:
% The mesh has NI cells in the i direction and NJ cells in the j direction.
% The x and y coordinates of the lower left hand corner of cell (i,j) are
% held in arrays x and y respectively which are both NI+1 by NJ+1. In
% this way the 4 vertices of cell (i,j) are (x(i),y(j)), (x(i+1),y(j)),
% (x(i+1),y(j+1)) and (x(i),y(j+1)). solid is an NI by NJ array which
% flags solid cells. If cell (i,j) is solid then solid(i,j)=1 otherwise
% solid(i,j)=0.
%
[m,n]=size(x);
NI=m-1; % number of cells in i direction
NJ=1; % number of cells in j direction
%
% Plot the mesh
hold on % keeps all plots on the same axes
% draw lines in the i direction
for i=1:NI+1
plot(x(i,:),y(i,:))
end
% draw lines in the j direction
for j=1:NJ+1
plot(x(:,j),y(:,j))
end
title('computational mesh')
xlabel('x')
ylabel('y')
% Fill in solid cells
for i=1:NI
for j=1:NJ
if (solid(i,j)==1)
solidx=[x(i,j),x(i+1,j),x(i+1,j+1),x(i,j+1),x(i,j)];
solidy=[y(i,j),y(i+1,j),y(i+1,j+1),y(i,j+1),y(i,j)];
fill(solidx,solidy,'k')
end
end % of j loop
end % of i loop
end % fdrawmesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fplotresults(xcen,u0,uinitial,uexact)
% Verification: not verified
% Dispays results on regular Cartesian mesh
[NI,NJ]=size(xcen);
u0=u0(2:NI+1,1); % extract computational values
plot(xcen,uinitial,'k--',xcen,uexact,xcen,u0,'k.')
%plot(xcen,u0)
title('plot of initial profile (--) and numerical(.) and exact solutions (solid line)')
xlabel('x [m]')
ylabel('U [kg/m]')
end % fplotresults
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%